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Abstract 



We present results from an extensive analytic and numerical study of a two- 
dimensional model of a square array of ultrasmall Josephson junctions. We 
include the ultrasmall self and mutual capacitances of the junctions, for the 
same parameter ranges as those produced in the experiments. The model 
Hamiltonian studied includes the Josephson, as well as the charging, Ec-, 
energies between superconducting islands. The corresponding quantum par- 
tition function is expressed in different calculationally convenient ways within 
its path-integral representation. The phase diagram is analytically studied 
using a WKB renormalization group (WKB-RG) plus a self-consistent har- 
monic approximation (SCHA) analysis, together with non-perturbative quan- 
tum Monte Carlo simulations. Most of the results presented here pertain to 
the superconductor to normal (S-N) region, although some results for the 
insulating to normal (I-N) region are also included. We find very good agree- 
ment between the WKB-RG and QMC results when compared to the exper- 
imental data. To fit the data, we only used the experimentally determined 
capacitances as fitting parameters. The WKB-RG analysis in the S-N region 
predicts a low temperature instability i.e. a Quantum Induced Transition 
(QUIT). We carefully analyze the possible existence of the QUIT via the 
QMC simulations and carry out a finite size analysis of Tqujt as a function 
of the magnitude of imaginary time axis L^-- We find that for some relatively 
large values of a = ^ (1 < a < 2.25), the Lr ^ oo limit does appear to 
give a non-zero Tqjjit, while for a > 2.5, Tquit = 0. We use the SCHA 
to analytically understand the Lj- dependence of the QMC results with good 
agreement between them. Finally, we also carried out a WKB-RG analysis 
in the I-N region and found no evidence of a low temperature QUIT, up to 
lowest order in a~^. 

PACS numbers: 74.20.-z, 74.50. -Fr, 74.60.Bt, 74.60.Dw 
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I. INTRODUCTION 



The physics of Josephson junctions arrays (JJA) has been a subject of significant interest 
in the last ten years . A large number of studies, both experimental 0-0 and theoretical 
||8|^3| have been devoted to them. Initially, part of the interest in JJA came from their close 



relation to one of the most extensively studied theoretical spin models, i.e. the classical 2-D 
XY model for which JJA give a concrete experimental realization. 

In non-dissipative JJA the two main contributions to the energy are the Josephson cou- 
phng between superconducting islands due to Cooper pair tunneling, and the electrostatic 
energy arising from local deviations from charge neutrality. In the initial experimental stud- 
ies, the size of the islands was large enough so that the charging energy contributions were 
very small, thus making the arrays' behavior effectively semi-classical. Recent advances in 
submicron technology have made it possible to fabricate relatively large arrays of ultrasmall 
superconducting islands separated by insulating barriers. These islands can have areas of 
the order of a few fim^, with self capacitances Cs ~ 3 x 10~^ fF, and nearest neighbors' 
mutual capacitance ~ 1 fF 0- Note that the mutual capacitance can be at least two 
orders of magnitude larger than the self capacitance. Therefore, these arrays have charging 
energy contributions, Ec, large enough so that quantum fluctuation effects are of paramount 
importance. 

In the Delft and Harvard ^ experiments, the island sizes were kept constant, while 
varying the normal state junction resistance, which in turn changes the Josephson coupling 
energy, Ej. This allows one to obtain arrays with values of the quantum parameter 

«m = (1) 

in the range [0.13-4.55] 0, or values as high as 33 [|]. In this equation we have used the 
definition of charging energy, 

ECra = (2) 

The experimental systems can be modeled by a quantum generalization of the classical XY 
model, because the phase of the order parameter associated with each one of the islands is 
canonically conjugate to its excess Cooper pair number. The magnitude of am determines 
the relevance of the quantum fluctuations. For small the quantum fluctuations of the 
phases are small and the system is well modeled by a renormalized classical 2-D XY model. 

The nature of the phase transition in the classical 2-D XY model is well understood, 
whereas in its quantum mechanical generalization there still are unsettled issues. One of 
the most notorious of these is the possibility of having a low temperature instability of 
the superconducting state. A possible reentrant transition was originally found within a 



mean field theory treatment of the self-capacitive XY model ||T6|J20| --|22| . An explicit two- 
dimensional study of the self-capacitative XY model, within a WKB renormalization group 
(WKB-RG) analysis also found evidence of a low temperature reentrant instability, triggered 
by a quantum fluctuation induced proliferation of vortices [T^ . 



Recently, Kim and Choi have studied the quantum induced fluctuations in these arrays. 



using a variational method p3[. They found that there is a range of values of the ratio 
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of charging to Josephson energy, for which there is a low temperature reentrance from a 
superconducting to a normal state. Similar results had been obtained by Simanek, also 
using a variational calculation, see for example Ref. ||22|| . 

A non-perturbative quantum Monte Carlo study of the self-capacitive model found a low 
temperature transition, but between two superconducting states The fully frustrated 
version of this model was also studied by quantum QMC and it yielded a larger jump 
discontinuity in the superfiuid density as compared to the one in zero field as well as the 
critical temperature one order of magnitude higher . A more recent analysis of the WKB- 
RG analysis has shown that, to lowest order in the quantum fluctuations, it must have the 
same critical temperature for a quantum induced phase transition (QUIT) A recent 

QMC study of the fully frustrated self-capacitive model by Mikalopas et al. |^ has suggested 
that the unusually large jump in the superfiuid density is dominated by metastability effects 
due to the particular nature of the excitations in the frustrated model. This result is in 
agreement then with the reanalysis of the RG equations. However, this study was carried 
out at relatively high temperatures and the question about the existence of a QUIT, both 
in the frustrated and unfrustrated cases remains open. We deal extensively with the later 
question here. Other studies find within MFT that to have reentrance it is necessary to 
include off-diagonal capacitances |[T^, while others do not agree with this finding PJT7|-pi? 



The search for a reentrant type transition is encouraged by some evidence of low temperature 
instabilities found experimentally in arrays of Josephson junctions M, ultrathin amorphous 



films ||2^, a multiphase high-Tc system [^], and in granular superconductors 

Most theoretical studies have been carried out using the self-capacitive model and dif- 
ferent kinds of MFT or self consistent harmonic approximations (SCHA) |19[. As already 
mentioned, these studies do not agree among each other on some of the properties of the 
phase diagram, in particular about the possible existence of a low temperature instability 
of the superconducting state. No study has been carried out that closely represents the 
experimental systems where both the self and mutual capacitances are explicitly included. 
The goal of this paper is to consider a model that is expected to represent the characteristics 
of the Delft experiments. In particular, we concentrate on calculating the phase diagram 
using different theoretical tools. 

One of the main results of this paper is presented in Fig. [l| which shows the vs. 
T phase diagram for an array with Cm > Cg both for the unfrustrated (/ = 0) and fully 
frustrated (/ = 1/2) cases. The left hand side of this diagram shows the superconducting 
to normal phase boundary (S-N) as data points with error bars joined by a continuous line. 
These data points were calculated using a QMC method, to be described later in the paper. 
We also show (as squares) the experimental results taken from Refs. For / = at 

small values of a^, the theoretical and experimental results agree quantitatively quite well 
with each other and with the semiclassical WKB-RG approximation. On the other hand, 
they only qualitatively agree for the / = 1/2 case and on the superconducting to insulating 
phase boundary. The normal to insulating transition line is shown to the right of the phase 
diagram. The latter is just a tentative boundary since our numerical calculations were not 
reliable enough to give the definitive location of this line, as also happens in experiments. 
The error bars in the calculated points used to draw the N-I line represent a crude estimation 
of the region where the inverse dielectric constant is different from zero. However, the issue 
of convergence of the calculation to the path integral limit is not resolved by these error 
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bars. As we will explain in the main body of the paper we found further evidence for a 
low temperature instability of the superconducting state in our numerical calculations. We 
found that this instability depends strongly on the magnitude of ctm and the finite size of 
the imaginary time axis in the QMC calculations. The latter finding sets strict constraints 
on some of the reentrant type behavior found in previous theoretical studies. Other studies 



have found reentrance very close to the superconducting to insulating transition |2^. This 
possibility is harder to study from our Monte Carlo calculations. 

The physical content of the phase diagram is generally understood in terms of the in- 
terplay between the Josephson and charging energies. For small a-^ and high temperature 
the spectrum of excitations is dominated by thermally excited vortices, which drive the 
superconducting to normal transition as the temperature increases, while the charging en- 
ergy contributes with weak quantum fluctuations of the phases. The latter produces, after 
averaging over the quantum fluctuations, an effective classical action with a renormalized 
Josephson coupling that lowers the critical temperature [11^ . 



For large and low temperatures the charging energy dominates. The excitations in 
this limit are due to the thermally assisted Cooper pair tunneling that produces charged 
polarized islands. At low temperatures, there is not enough thermal energy to overcome 
the electrostatic coulomb blockade so that the Cooper pairs are localized and the array is 
insulating. As the temperature increases, the electric dipole excitations can unbind, driving 
an insulating to conducting transition (I-C). In the hmit ^ Cs, it was suggested that 
the I-C transition could be of the Berezinskii-Kosterlitz-Thouless (BKT) type, for in this 
case the interaction between charges is essentially logarithmic p9| , p0| ,p|. However, for the 
experimental samples it has been shown that rather than a true phase transition what is 
measured appears to be a crossover between an insulating to conducting phase, characterized 
by thermally activated processes [^,§1. It is likely that the reason for the crossover is the 
short screening length present in the samples (A ^ 20 lattice sites). Both experimental 
groups find that a simple energetic argument gives an explanation for the activation 
energy found in the experiments. Furthermore, the nature of this crossover may be linked 
to thermal as well as to dynamical effects. As we shall see, theoretically the I-N phase is 
hard to study in detail. 

For large values of the model can be approximated by a 2-D lattice Coulomb gas, 
where a perturbative expansion can be carried out using Ej as a small parameter. This type 
of calculation was performed in Ref. |Q in the limit Cs -C Cm. The analysis lead to a 2-D 
Coulomb gas with a renormalized coupling constant. Here we will extend this calculation 
to obtain a more accurate estimation of the renormalized coupling constant. We do this 
because we are interested in seeing if it is possible to have a QUIT instability in the low 
temperature insulating phase. This possibility is suggested by the dual symmetry of the 
effective action between charges and vortices found in Ref. and the fact that the am 
perturbative expansion shows a low temperature QUIT instability in the superconducting 
phase. We find that the results of a first order expansion in a"^ does not present this type 
of low temperature instability. 

Among the most interesting regions of the phase diagram is when the Josephson and 
charging energies are comparable. For this nonperturbative case, using a path integral 
formulation and the Villain approximation [^], an effective action for logarithmically inter- 
acting charges and vortices was derived in Ref. P in the case where Cg -C Cm. The action 
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of the two Coulomb gases shows an almost dual symmetry, so that at an intermediate value 
of Om, both the S-N and the I-C transitions converge to a single point as T — 0. A similar 
picture was derived in Ref. using a short range electrostatic interaction and a mean field 
renormalization group calculation. A nonperturbative calculation is needed to determine 
the actual shape of the phase diagram in this region. We further discuss this point in the 
main body of the paper. 

It has been argued that at T = the self-capacitive model is in the same universality 
class as the 3-D XY model where the ratio as = {q^ /2Cs) / Ej would play the role 

of temperature. This analogy would result in a transition at some finite value of from 
a superconducting to an insulating phase. Moreover, there should be a marked signature 
in the nature of the correlation functions when crossing over from a 2-D XY model at high 
temperatures to a 3-D XY model as T — 0. When as = 0, the correlation functions decay 
algebraically in the critical region, with a temperature dependent exponent. At T=0 and 
as ^ there is a single critical point at as = ac, so that the correlations decay exponentially 
above and below ac, and algebraically at as = ac- The question is then, how do we go from 
algebraic to exponential correlations as T changes? This can only happen by having a change 
of analyticity in the correlations, thus the possibility of having a QUIT in the self-capacitive 
model. 

The situation is different in the mutual-capacitance dominated limit, of experimental 
interest. When Cs = the model is equivalent to having two interacting lattice Coulomb 
lattice gas models. The general critical properties of this case are not fully understood 
at present. A further complication arises when Cg is small but non-zero. The map to a 
higher dimensional known model does not work in this case, and the problem has to be 
studied on its own right. Because of the essential differences between the self-capacitive 
and the mutual-capacitance dominated models one can not just take results from one case 
and apply them to the other. It is the goal of this paper to explicitly study the mutual 
capacitance dominated model, but with nonzero Cg- A brief report on some of the results 



of this paper has appeared elsewhere [24 



The outline of the rest of the paper is the following: In Section |l| we define the model 
studied and derive the path integral formulation of the partition function used in our calcu- 
lations. In Section |III A| we present a derivation of the path integral used in the semiclassical 
analysis, in the limit where the Josephson energy dominates. We carry out a WKB expan- 
sion up to first order in a^, finding an effective classical action where the charging energy 
contributions are taken into account as a renormalization of the Josephson coupling. In 
section [IIIB| we find general renormalization group (RG) equations from which we obtain 
the phase diagram for small Om- In Section |V] we study the large a^ limit, in which the 
charging energy is dominant over the Josephson energy. There we obtain an effective 2-D 
Coulomb gas model with a quantum fluctuations renormalized coupling constant. In Section 
V A| we discuss our QMC calculations and define the physical quantities calculated. In Sec- 
tion |VB| we present some technical details of the implementation of the QMC simulations. 
In Section |V C| we give the QMC results for / = mainly, but also for / = 1/2. There we 
make a direct comparison between the semiclassical approximation results, the QMC calcu- 
lations and experiment which lead to the phase diagram discussed above. In that section 
we also present an Lr dependent analysis of the apparent Tqujt for three relatively large 
values of am = 2.0, 2.25 and 2.5. The ^ oo extrapolation of the results leads to a finite 
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Tquit for am = 2.0 and 2.25, while for am = 2.5 we get a Tquit{Lt- = oo) = 0. In section 
[V^ we discuss a self-consistent harmonic approximation (SCHA) analysis, that we use to 
analytically study the phase diagram, and that helps us understand the finite size effects of 
the imaginary-time lattices studied in the QMC calculations. At the end of the paper there 
are two appendices where we give more technical details of the analysis. In Section |V1I| we 
restate the main results of this paper. 



II. THE MODEL AND THE PATH INTEGRAL FORMALISM 

In this section we define the Josephson junction array model considered in this paper 
together with the path integral formulation of its corresponding partition function. 

We assume that each superconducting island in a junction can be characterized by a 
Ginzburg-Landau order parameter ^'(r) = \'^Q(r)\e'^'^^^^ , where r is a two-dimensional vector 
denoting the position of each island. If the coherence length of the Cooper pairs is larger 
than the size of the islands, we can assume that the phase of the order parameter is constant 
in each island. Moreover, the amplitude of the order parameter is expected to have small 
fiuctuations about an electrically neutral island and can then be taken as constant trough 
the array. We will assume that the charge fiuctuations have an effect on the electrostatic 
energy but not on the Josephson contribution to the Hamiltonian. The gauge invariant 
Hamiltonian studied here is 

n = Hc + Hj = ^ J2 niri)C~\n,f2)nif2) + Ej J2 [l - cos (0(fi) - ^(rl) - A,-,,,-,)] , 

<fi,r2> <fi,f2> 

(3) 

where q = 2e; 0(r) is the quantum phase operator and n(r) is its canonically conjugate 
number operator, which measures the excess number of Cooper pairs in the r island. These 
operators satisfy the commutation relations [n{ri) , (f){r2)] = ~iSfi,f2 HH- Here A^i,^*^ 
defined by the line integral that joins the sites located at ri and r2, Af^^^^f^^ = ^ JI^ A-dl, where 

A is the vector potential and $o is the fiux quantum. In Eq. He is the charging energy 
due to the electrostatic interaction between the excess Cooper pairs in the islands. The 
C~^(ri,r2) matrix is the electric field propagator and its inverse, C(ri,r2), is the geometric 
capacitance matrix, which must be calculated by solving Poisson's equation subject to the 
appropriate boundary conditions. This is not easy to do in general and typically this matrix 
is approximated, both theoretically and in the experimental analysis of the data, by diagonal 
plus nearest neighbor contributions [^: 

C(n, r,) = (a + zC^)5,,,,, -C^Y. K^ur^+d- (4) 

d 

Here the vector d runs over nearest neighboring islands, z is the coordination number, Cg 
is the self-capacitance of each island, and Cm is the mutual capacitance between nearest 
neighbor islands. In the experimental arrays, typically Cm ~ lO^Cg ~ IfF [Q. 

The second term in Eq.@ is the Josephson energy, which represents the probability of 
Cooper pair tunneling between nearest neighboring islands. The Josephson coupling energy 
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Ej = $o^c/(27r) is assumed to be temperature independent, where ic is the junction critical 
current and $o the flux quantum. Here we are interested in calculating the thermodynamic 
properties of the model defined by 7i. The quantity of interest is the partition function 



Z^Tr{e-'3^}. 



(5) 



The trace is taken either over the phase variables, 0, or the numbers operator, n. To evaluate 
the partition function we will use its path integral representation To derive the path 

integral we use the states 



< n(fi)|0(f2) >= (5^-1 



exp{in{fi)(f){fi)} 



)"2 



(6) 



We will also use the fact that both {|n(f^ >} and {\(p{r) >} form complete sets. To start 
we write the partition function as a trace in the phase representation 



Z = Y[J^^d(l){0,r) < {0(O,r)}|exp{-/37i}|{0(O,f)} 



> . 



(7) 



As usual we use the Trotter formula 
exp{-f3{Hc{n)+Hj{4>))} = 



Lr 



ex^{-{(3/Lr)Hc{h)]ex^{-{(3lLr)Hj{<P)} 
+0 (l/L^) . (8) 

Next we introduce Lr — 1 complete sets {|0(t, r) >}, r = 1, 2, . . . , L,- — 1 in Eq. (|^) so that 



Z = nn / t^0(r,r)<{0(O,r-)}|exp{-(/?/L,)7^}|{0(l,f)}> X 

r T=0 

X <{0(l,f)}|exp{-(/5/L,)7i}|{0(2,r)}> X 

{<\>{U - 1, r)}| exp {-{(5lU)n] |{0(O, r)} 



X ■ ■ ■ X < 



> 



+0 

At this point we need to calculate the short time propagator. 



(9) 



< {</)(r, r)}|e-(^/^-)^|{0(r + 1, f)} >= E < r)}|e-('^/^-)^|{n(r, f)} > x 

n(T,r)=— oo 

X <{n(r,f)}|{0(r + l,f)}>, (10) 

where we used a summation over the complete set |{?t-(t, r)} >. From Eqs. (^) and (H) this 
propagator can be written as 



1 



< {0(r,f)}|e-('^/^^)^|{0(r + l,f)} >=Viir- E ^ 



i n(r,r)[(/)(T+l,r)— ()!>(r,r)] 



r n(T, rj=—oo 



g-(/3/L0W({n(r,r)},{<^(r,r)}) ^ 



+0(1/L^). 



(11) 
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Inserting this equation in Eq. (j^) we obtain the following path integral representation of 
the partition function 



^= nn 

T=0 f 



X exp 
+0(1/L^). 



j-2TT 


d(f){T, 


Jo 


271 


Lr 


E 

T = 



E exp 

{n{r,r)}=— oo 



J2 n{T, f) [0(r + 1, f) - 0(r, f)] 



r=0 

— n 
2 



+ 



(12) 



together with the important boundary condition (f){Lr, r) = 0(0, r). These equations are our 
starting point for the semiclassical approximation analysis discussed in the next section. 



III. WKB AND RENORMALIZATION GROUP EQUATIONS 

A. Semiclassical limit 

The semiclassical limit corresponds to taking 0, or 0. The summations 

over {n(r, r)} in Eq. ( [T^ ) can be carried out and the result leads to 0(r + l,r) = 0(r, r), 
for r = 0, 1, . . . , Lt- — 1. In other words, in this limit all the phase variables are constant 
along the imaginary time axis, and we recover the classical 2-D XY model [^,^. As the 
charging energy increases the value of 0(r, r) fluctuates along the r-axis; these fluctuations 
suppress the XY phase coherence in the model lowering its critical temperature. For the 
self-capacitive model (Cm = 0), at T = 0, one can map the model to an anisotropic three- 
dimensional XY model H,pJ^. This model should have a transition between ordered and 
disorder phases at a critical coupling {EcjEj)^. Here Ec^ = e^/2Cs, so we would expect 
the phase boundary to go all the way down to T = for large enough charging energy. 

In this section we study the change of the critical temperature as Ec^ increases, for small 
values of the ratio = Ec^_^J Ej. We start by eliminating the {n's} the from Eq.(0) using 
the Poisson summation formula 



E /H = E 

n=— oo 



m=—oo 



/(x)e'~da; 



obtaining 



Z=nv^n/„ J^^Kr) E exp 



r=0 



Here we defined the action 



27rpq^ 



{m(r,r)}=— oo 



SmAm}] 



(13) 



(14) 



Lr-l 



smAm}]= E 



T=0 L 



ri,r2 



xC(fi,f2)[0(r+l,f2) 
+0{1/Ll). 



(r, r2) + 27rm(r, f2)] 



+ 



(15) 
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It is convenient to write the paths in the partition function separated into a constant part, 
that corresponds to the classical model, plus a quantum fluctuating contribution, over which 
we will perform the integrations to find an effective classical action. First we eliminate the 
summations over the {m's}. This is done at the same time that the integrals over the {4>'s} 
are extended from [0, 27r) to (—00,00). After a couple of standard variable changes p2|j3^ 
we get an action where the phases and the charges are separated. 



V m(fi)C(fi, f2)m(f2) + - / dr 

fi Jo 

ri,r2 



X 



X 



2g2 ^ dr 

^ ri,r2 



dr 



(16) 



Here the variables ip{(5h, r) = ilj{0,r), and the integers m{r) are called the winding numbers. 
This equation shows that the winding numbers are the charge degrees of freedom and that 
the coupling between phases and charges appears only in the Josephson term. We can also 
see from this equation that in the semi-classical limit (small charging energy) the charge 
fluctuations are exponentially suppressed. This is more so for the m's because they have 
a discrete excitation spectrum. Therefore, to lowest order in the semiclassical analysis we 
will set m(f^ = 0, leaving integrals only over the phases. Next we separate the ip's into a 
constant plus a fluctuating part 



-i/,(r, r) = (f){r) + </)/(r,r). 



(17) 

At this point we use the following argument ||35|,^. First that the Lagrangian is invariant 
under the transformation ip{0, r) —>■ tp{0, r)+27Tl{r) for all integers /(r), so that we can extend 
the limits of integration over ilj{0,r) to (—00,00), safe for an extra overall multiplicative 
constant. Now, the limits of integration for 0(r)e(— 00, 00), and because of the periodicity 
of (pf, we can Fourier series expand it as 



<Pf{T,r) = (/5;i)-i/2^[0fc(r)e*-'=- + aa], (18) 

k=l 

where the uJk = 2T{k/(3h are the Bose-Matsubara frequencies. We have then the partition 



function 33 



-. J — ( 



(27r/5g2)i/2 n 



ujlh 



/oo /"OO 
dRe0fc(f) / rflm0fc(r) 
-00 J—oo 



xexp^--5[{0},{0^}] 



(19) 



Next, we expand the Josephson term in the action up to second order in (j)f, for higher 
order terms are suppressed in the integrations. After performing the integrations over the 
Euclidean time r, and the Gaussian integrations, the effective partition function reads 
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(2vr/5g 

CX) 

k=l 



I ^ 

det <^ 5r^,,r'2 + T2-^ C"^(ri, ' 



(20) 



Here we want to expand this partition function in powers of the charging energy. This is 
equivalent to expanding in powers of g^, so our next step is to expand the determinant. We 
use the following identities 



det [I + D] = exp{Tr[ln(I + D)]}, 
ln(I + D) = -Y.- ^D", 



n=l 



n 



(21) 
(22) 



where I is the identity matrix. To lowest order in q and using the resuh E^i[gV(^^fc)^] 
(g/3)^/24, the effective partition function for (j) is then 



d(j){f) 



^ ^-U^ dHj 



24 



90 (fi) 90 (fa) 



(23) 



To further advance the calculation we now use the properties of the Josephson energy. 
We start by using the fact that it is a local nearest neighbor interaction, so from Eq. (|^) 



(24) 



d 



with the d running over the nearest neighbors to fin the lattice. From this equation we can 
see that the second derivative of ifj({0}) is given by 



90 (fi) 90 (fa) 



E 



r (0(fo - 0(fi - d)) {5,,,., - 5. . ^,-) + 

+r (0(ri + d) - 0(fi)) (5.1,^, - b,,^lr^ 



(25) 



where f"[x) = d'^f{x)/dx^. 

In this paper we consider a periodic array, which implies that the inverse capacitance 
matrix is invariant under translations and rotations, that is C~^(fi,f2) = C~^(|fi — fal). In 
particular, this makes C~^(f, f ± d) = C~^{\d\), independent of the direction of d. Notice 
that here we are using d to denote the vectors that connect nearest neighboring islands, 
therefore in a periodic and symmetric array all of them have the same magnitude, allowing 
us to take the terms C~^(|(i|) out of the summations. From these considerations the trace 
gives 
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ri,r2 



d(f){fi)d(j){f2) 



c-\\o\) - c-\\d\)] j: j: f"U{f+ d) - m) . (26) 



Next notice that since f"{x) = —f{x) (up to a constant), both terms in the argument of the 
exponential in Eq. ( p3D are the same cosine function of the classical phase variables, with 
only different coupling constants. Finally, the effective semi-classical partition function can 
be written as 



Zeff = n ydet[C] J 



exp{-/?,ff/f^({0})}. 



-oo (27r/3g2)i/- 



where the effective temperature is explicitly given by 

= f3 - q'^[C-\\% - C-\\d\) 



(27) 



(28) 



Notice that to obtain this result we have used an argument that could be questionable, 
namely the extension of the 0(0, r) to the (— cxo, oo) range in the path integrals. All the other 
approximations are consistent with the semiclassical approximation and the symmetries used 
are exact. However, to continue we will now restore the [0, 27r) range of the phases to use the 
results known from the BKT theory. As we will show later in the paper, the nonperturbative 
QMC results do agree quantitatively with the WKB results and experimental results to be 
discussed later. A similar effective result was first obtained but for the self-capacitive model 
in ig. 

One of the important properties of Eqs.(^) and ( p8l) is that up to this point we have made 
no assumptions about the structure of the capacitance matrix that go beyond translational 
invariance. Later on we will make specific choices of this matrix when we make direct contact 
with experimental findings H]. 



B. Renormalization group analysis 

Now that we have expressed the quantum mechanical problem as a modified 2-D classical 
XY model we can directly apply the well known results for this model 0,^]. The standard 
physical picture of the excitation spectrum in this model is of spin-waves plus vortex pair ex- 
citations. At low temperatures the energy to create an isolated vortex grows logarithmically 
with the size of the system, therefore excitations are created as bounded vortex-antivortex 
pairs. As the temperature increases, the vortex pair density increases until they unbind 
at a critical dimensionless temperature Tbkt = 0.894(5) []37| , ^ . The BKT scenario is best 
understood in terms of a renormalization group (RG) analysis [ P^p5| . The RG flow diagram 



is obtained from a perturbative expansion in powers of the vortex pair fugacity y. To lowest 
order in the RG equations corresponding to our problem are 

^ = -A^KW. (29) 
^ = [2 - 7rKeff]y. (30) 
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Here we have used the following definitions: 



Kes = K- xK' 

-,2 



X 



12Ej 
K = (3Ej. 



C'\\0\) - c-\\d\) 



Then the equations for the coupling constants K and y are 

dl ^ ^ {2Kx-iy 

^ = [2 - nK{l - xK)] y. 



(31) 
(32) 
(33) 

(34) 
(35) 



To find the critical temperature, we use the initial conditions from the temperature and the 
bare vortex pair fugacity 



K,s{l = 0) = l3,^Ej 
y{l = 0) = exp <^ - 



1 + 



1 



O 



1 



TC 



K,s{l = 0) 



(36) 
(37) 



The RG equations have two nontrivial fixed points (for x < vr/8). One corresponds to the 
effective BKT thermal fluctuations driven transition, and the other to a quantum fluctuations 



induced transition (QUIT) [|T3]-[T5|. 

One way to analyze the structure of the RG flow in the {y, K) phase space is to use a 
conserved quantity associated with Eqs. (|29| ) and ( pO] ) 



A = -vrlnK 



eff 



K. 



+ 2TT^y^. 



eff 



Using Eq. (|32D and expanding up to flrst order in x we flnd 



A = TixK 



TflnK — — 
K 



2vry. 



(38) 



(39) 



Figure |^ shows the RG flows obtained from numerically solving the RG equations for 
different values of A, where the arrows indicate the direction of increasing /. We have also 
plotted the set of initial conditions from Eqs. (^6]) and (^) as a discontinuous fine. One 
important flow line is the separatrix between the lines for which y{l ^oo) — >■ oo and those 
for y{l—>-oo) 0. This line has Ac = — 7r[l + ln(2/7r)], which is determined by the condition 
that it must pass through the point {y = 0, K^s = 2/7r). The critical temperature is obtained 
from the intersection of the separatrix with the initial conditions given in Eqs. (|36|) and 
(p7|). This intersection exists only if x is less than the critical value Xc < tt/S, to be estimated 
below. Fortunately, we do not need to flnd this intersection explicitly since we already know 
the critical value of the effective coupling K^s, which is the usual critical coupling of the 
classical XY model, K^^ = K^^ ^ 1.1186 Therefore, the values of the two critical 

couplings are 
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K{l-xK) = Kf\ (40) 
1 



K± = — 
2x 



1 ± a/1 - AxK^^ 



(41) 



which leads to = l/(4irf^) ^ 0.2235. 

We note from Fig. 2 that the axis can be divided into three different regions. If 
we set = Kqun^, and K_ = K^kt, then in the region [-R'quit' -^bkt]; ^ increases, the 
fugacity of the vortex-antivortex pairs decreases. In the limit / oo the energy to create 
a macroscopic vortex pair becomes infinite. Therefore, the system is superconducting for 
temperatures in this interval. For temperatures T > EjK^]^j. the renormalized fugacity 
y{l) increases and the low y approximation breaks down. For these temperatures, the array 
is normal. For T < EjKqjjjrp, the vortex pair density increases due to the quantum fluc- 
tuations, leading us to think that there may be a low temperature transition driven by the 
quantum fluctuations (QUIT). 

To obtain the results described above we used a high temperature perturbative calcu- 
lation. Therefore the QUIT results are in principle outside the regimen of validity of the 
WKB-RG calculation. We need, then, other calculations and approaches valid at low tem- 
peratures to prove or disprove the existence of the QUIT. Expanding Eq. (^), up to first 
order in x and using Eq. (|33|) we find the critical temperatures. 



'-bkt'-^T^bIt-t^^ + 0{x'), (42) 

Tquit-'^x + 0{x^). (43) 
Kb 

Note that these equations are applicable not only in 2-D, for if the system described by 
Eq. (^) has a transition point at some K^^ then the equation K'^^^ = K — xK^ has two 
solutions for K. In this argument we should notice that the existence of the second solution 
for K depends on the higher order terms in the x expansion. Note that the change of Tbkt 
for small x is correctly given by the small x result. The existence of a low temperature 
quantum phase is of a nonperturbative nature, however. That is one of the reasons why we 
resort to using the nonperturbative quantum QMC approach latter in the paper. 

Another interesting property of Eqs. (H2|) and (^) is that the first order correction 
does not depend on the specific value of Tj^kt- particular, if we add a magnetic field 



to the Hamiltonian in Eq. (y), all the calculations leading to Eqs. (|27|) and (|28|) would be 
unchanged. Therefore, if Tj^^^B) is the superconducting to normal transition temperature 
for the array in a finite magnetic field i? at x = 0, then to first order in x we must have 

UB) ^ T^'\B) - {Ej/k^)x + 0{x^). (44) 



This equation is in agreement with the results obtained in Ref. |]I4[. Furthermore, we notice 
that to lowest order in x, the Tqujt must be the same with then without a magnetic field. 
This result was not noted before and it can be used as a test of the QMC calculations, in 



particular those of Mikalopas et al. |25]. 



To make comparisons with experiment, we need to specify the capacitance matrix. In 
particular, if we use Eq. (^ and the specific geometry of the array, we find that x is given 
by 
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X 



l2zEjC„ 
UzEjC, 



i-ac-i(lo|) 



(45) 



The function g{w) can be written as an elliptic integral for a two-dimensional square lattice 
|38| . For a general lattice geometry, we find the following limiting behavior 



z — z{\ + if w <^ 1, 



w -"^ {1 — (47rw) ^Inw}, if w ^ 1. 



(46) 



Using Eqs. 



I) and (^61) we get 



i^B-L BKT 



(2/3)as + 0(a2), if ^ > C^, 
(2/3;2)a^ + 0(O, if ^ « C^m- 



(47) 



Ej Ej 

This result is in agreement with the Monte Carlo calculation of the superconducting to 
normal transition temperature carried out in Ref. for the self-capacitive model. As we 
will show later in this paper, it is also in good agreement with our QMC calculations for the 
model dominated by the mutual capacitances. 



IV. INSULATING TO NORMAL CROSS OVER. 

So far we have studied the normal to superconducting transition in the limit where 
the Josephson energy dominates over the charging energy. In the opposite limit, when the 
relevant excitations are charge fluctuations, the transition is expected to be from a normal 
conducting state, where the charges are free to move, to an insulating state where the charges 
are bound into neutral dipole pairs (see Fig. |l]). It has been suggested that in the limit 
(Cs/Cra) — > this I-N transition would be of a BKT type P^ , P(]| flJ5[1 . Experimental results 
have shown, however, that the behavior of the fabricated samples is better explained by a 
crossover from a normal to an insulating phase [§H^ . In finite systems, like the experimental 
ones, we would expect a rounding of the transition. Furthermore, in the experiments, the 
screening length is shorter than the sample size A ~ \fc^JC's ~ 18 lattice spacings 0. 
Minnhagen et al. have argued that for any finite screening length, the transition is washed 
out even for an infinite array This is not difficult to understand since in the BKT 

scenario the superconducting to normal transition depends on the unscreened nature of the 
vortex logarithmic interaction p^|35 . 



In this section we present results from a perturbative calculation of the effect of the 
Josephson energy on the expected I-N crossover temperature. We start with Eqs. (p!6| ) and 
(p!7|) leading to 



r ' ' ^ m(rl=— oo r=i r 



m{f)=—oo ' 



;d(j)f{T,r) X 



X exp 



1 f^fi 



drLi 



(48) 
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where the action is now given by 



1 rP^ 



Jo 2 Dq^^ h Jo 



1 /-/^^ 



h Jo 



2 Pq 



X 



ri,r2 



2g2 dr 



dr 



+Hj{{(P{f) + (Pf{T, r) + (27r//5^)m(r)r}) 



with the boundary condition 



(49) 



(50) 



Since we are interested in the charge degrees of freedom, our task here is to integrate out 
the phases. This hmit has been studied before, in particular in Ref. §]. Here we are not 
only interested in the crossover temperature, but we mostly want to ascertain if there is an 
equivalent QUIT in the insulating phase at low temperatures. 

Since we are at the limit Ej <^ Ec, the Josephson energy can be treated as perturbation. 
We expand the exponential 



exp 



1 rP^ 



h Jo 

We note that Eq, 



drHjir) 



1 - - / drHjir) + / dTdT'Hj{T)Hj{r') + ... (51) 

n Jo 2h Jo Jo 



I i-i3h i-pn 



f m(r)=— oo 



can be written as 
(2vr)2 



2/3g2 



m(ri)C(fi, f2)m(r2) 



ri,r2 



Zcs{{m}), 



(52) 



where Z^ contains only phase degrees of freedom and can formally be written as 

1 



/oo 
V(j)f{r)exp 
-oo 



h 



-s 



2q' Jo 



I3h 



d(f)f 



d(f)f 



ri,r2 



(53) 
(54) 



Here we have used the following short hand notation for the measure 



Lt-1 



P0^(f)= hm nVdet[C]n 

T=L r 



Lr 



;d(t)f{ 



(55) 



noting that strictly speaking the integrals over a finite number of L^-'s have to be calculated 
before the limit L,- — > cxd is taken [^ . 

All the interactions between phases and charges are contained in Zeff({m}), the effective 
partition function for the charges. The details of the explicit evaluations of Zf,^ are given in 
Appendix A. The result for Eq. (|5^ can then be written, up to second order in Ej, as 
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r m(r}=— oo 



1 



1 /a 



(m(fi) — mfro))^ — , ^ 



+ 



2 \ g2 



x(m(fi)-m(f2),K[l-aC-i(|0| 

/ <ri,r2> 



+ 



Here we have defined 

k = 

J(m, k) = 



(27r)2a 

1/2 



(ia;i cos(27rmxi) exp — {2(27r) — xi) 



(56) 

(57) 
(58) 



The function T(m, K) is an even function of m, so it can be expanded in a Taylor series 
in m^. One way to do it is to take cos(a;i) ?a 1 + (l/2)x^ — (l/24)a;f + . . .. This is a good 
approximation if the coefficient in the exponential is large. Since we are interested in discrete 
values of m, and considering that values of m greater than one are suppressed even near the 
transition point, we can use the following approximation 



j(m, k) = j(o, k) - [j(o, k) - k) 

With this approximation we can write Eq. (E6|) as 



m 



f m(r)=— oo 



1 



E 



^ (m(ri)-m(f2))2--^ — ^^m(r) 

/Acff A^ 2a VC-n 



<ri,r2> 



The effective coupling constant is given by 



h{w) = 



1 + 

1/2 



h{K[l-aC-'{\0 



— cos(27rx)] exp — {2(2nY / z}w x{l — x) . 
The function h{w) has the following limiting asymptotic behavior 

{l/2)w^\l - w (1/-2)(12 + (27r)V3)l + 0{w'^), if < 1, 



(59) 

(60) 

(61) 
(62) 



h{w) 



+ j^{z/2)w~' - j^{zl2f w-^\ + 0{w-% if w » 1. 



(63) 



We now use the fact that for the experimental systems Cs ^ Cm, so that in the limit 
(Cs/Cm) — i> we can use 



lim 

(a/Cm)^o 



i-c^c-^doi) 



lim 

(Ce/CnO^O 



1 + 1_ In 

Ait \Cm/ VCm 



1. 



(64) 
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The end result is a discrete Gaussian model with an effective coupling constant given by 
Eq. (|6l[) . This effective model can be transformed into a Villain model |40,43|. The critical 



points of this model are given by the equation K^s = , where is the critical coupling 
for the Villain model, ~ 0.752(5) for a square array [Q. In other words we have to 
solve the equations, 

K, = Kj + nh{K,), (65) 
n = Kj (7rV4) {Ej/Ecy. (66) 

From these equations we find the first order correction to the crossover temperature for a 
square array 

'^^ - 0.259 f^Wo((Ej/Ec;)^). (67) 



Here we have used h{K^) ~ 0.0106. 

An important property of Eq. (|65|) is that we can show that it has only one solution, 
since the function h{K) is concave for small K and it has an inflection point at i^infi, 

A sufficient condition for Eq. (|65D to have only one solution is > Ki^^. This condition is 
satisfied for square as well as triangular arrays. This result shows that there is no insulating 
QUIT phase and it is in clear contrast to the existence of the QUIT found using the WKB- 
RG approximation in the superconducting phase. 



V. QUANTUM MONTE CARLO RESULTS 

A. Definition of Physical Quantities Calculated 

The two important physical parameters in our analysis are the temperature and a = 
Since in the experiments the self-capacitance is much smaller than the mutual capacitance, 
the relevant quantum parameter here is 

Er 

In the region where am is small, the phases dominate and we expect a superconducting to 
normal transition. The quantity we will use to characterize the coherent superconducting 
phase is the helicity modulus defined as 



T 



(70) 

A=0 



Here x is the unitary vector in the x direction. The superfluid density per unit mass, ps, is 
proportional T, with Ps{T) = y (^) '^{T), where a is the distance between superconduct- 
ing islands, m is the mass of the Cooper pairs, and V is the volume. From Eqs. (|12[) and 
m we get 
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EjLxLy 



-T(r) 



1 



LxLyLj- 



ILt-I 



Lr 



T=0 r 

/Lr-1 



T=0 T,f I 



(71) 



The quantity we shall use to probe the possible charge coherence in the array is the inverse 



dielectric constant of the gas of Cooper pairs, defined as ||47| , ^ 

g2 1 



- = lim 



< n{k)n{—k) > 



(72) 



ksT C{k) 

We can obtain the Fourier transform C{k) from Eq. (^) for the capacitance matrix to get, 
C{k) = Cs + 2Cm[l - cos{kx)] + 2Cm[l - cos(A;y)]. (73) 
The Fourier transform of the charge number is defined by 



n{k) 



LxLy r 



n(r) exp 



ik ■ r 



(74) 



Using this equation we can obtain a path integral representation for this correlation function, 
given by 



< n(ri)n{r2) >- 



^ r=0 r 



2tt 



{m(T,f)}=—oo 



exp 



d(j){Lr,fi)d(j){Lr,f2 

The action is given in Eq. (p!5|). This equation becomes 



Sm,{m}] 



(75) 



0{L,,r)=0(O,r) 



r=0 

1 dS 1 dS 



X 



^si{'#'}.Wl 



exp 



(76) 



with 



1 dS 



hd(f){Lr,fi) Pq^ 



J2 C(ri, f) [0(L„ f) - 0(L. - 1, f) + 27rm(L. - 1, f)] . (77) 
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Notice, that this is not a well behaved operator since in the limit L,- ^ oo we would have 
to subtract two large numbers and the path integral in the second term in Eq. (^) would 
diverge. This divergence is canceled out by the first term in Eq. ([761). This can be seen 
explicitly by doing the calculation of seting Ej = 0, which leads to 

< n{fi)n{f2) >= C(ri,f2) + {^^J~^ XI C(fi, r3)C(f2, r^) < m{f^)m{ri) > . (78) 

This result can be put into Eq. (|72D to obtain a finite inverse dielectric constant, 



1 

e 



lim 

fc-+0 



C{k) < \m{k)\' > 



(79) 



Here we have used the Fourier transform defined in Eq. ([74[) and the m(r) defined as 



T=0 



m(T,r] 



Note that in general this operator will not exactly be the inverse 
dielectric constant of a gas of Cooper pairs, since it will depend on L^. But we expect that 
it does contain most of the relevant information of the inverse dielectric constant of our 
charged system. In our Monte Carlo calculations we have used the general result Eq. (|76D 
valid for Ej ^ and finite L^-. 



B. The Simulation Approach 

Up to now we have seen that the partition function defined by the Hamiltonian in Eq. 
(§) can be expressed in different convenient representations for analytic analyses. To carry 
out our QMC calculations, we have used what is, in principle, the most straightforward 
representation of Z given by Eqs. (|14D and ([15[) ; it involves the phases and the charge 



integer as statistical variables. This representation is general enough to be used over all the 
whole parameter range covered in the phase diagram. 

In this case we have a set of angles 0(r, r) G [0,27r), located at the nodes of a three- 
dimensional lattice, with two space dimensions, and Ly, and one imaginary time dimen- 
sion, Lt-. The periodic boundary condition, comes from the trace condition in Eq. (^, and 
we also have chosen to use periodic boundary conditions in both space directions. The link 
variables m(r, r) are defined in the bonds between two nodes in the r direction and they 
can take any integer value. 

We have basically used the standard Metropolis algorithm to move about in phase space 
As the phases are updated we restrict their values to the interval [0, 27r). Moreover, the 



shifts along a r-column and the individual phase moves are adjusted to keep the acceptance 
rates in the range [0.2,0.3]. 

If am is small, the system is in the semiclassical limit. In this case the fluctuations of the 
phases along the imaginary time axis as well as the fluctuations in the m's are suppressed 
by the second term in Eq. (|1^) . Attempts to change a phase variable will have a very small 
success rate. Therefore we implemented two kinds of Monte Carlo moves in the phase degrees 
of freedom. In one sweep of the array we update the x Ly imaginary time columns, by 
shifting all the phases along a given column by the same angle. This move does not change 
the second term in Eq. (0), and thus it probes only the Josephson energy [[1^]. To account 
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for phase fluctuations along the imaginary time axis, which become more likely as (am/T) 
increases, we also make local updates of the phases along the planes. 

Another aspect of the implementation of the QMC algorithm is the order in which 
we visit the array. This is relevant for the optimization of the computer code in different 
computer architectures. In a scalar machine we have used an algorithm that updates column 
by column in the array. For a vector machine we have used the fact that for local updates, 
like the ones we use, the lattice can be separated into four sublattices in a checkerboard-like 
pattern. This separation is done in such a way that each of the sublattices can be updated 
using a long vector loop without problems of data dependency. Using this last visiting 
scheme, the cpu time grows sublinearly with the size of the array. One of the problems that 
this type of visiting scheme has in a vector machine, like the Cray C90, is that the array's 
dimensions have to be even, and this produces memory conflicts. We have not made attempts 
to optimize this part of the code. We have not used parallel machines in our calculations 
but the same type of checkerboard visiting scheme would lead to a fast algorithm. 

We followed Ref. [|I4[ and replaced the U(l) symmetry of the problem by a discrete Z^v 



subgroup. We took = 5000. This allows us to use integer arithmetic for the values of the 
phase variables, and to store lookup tables for the Josephson cosine part of the Boltzmann 
factors. This simplification can not be used for the charging energy part of the Boltzmann 
factors, except in the = case, where the m's can be summed up in a virtually exact 
form. In the latter case we can also store lookup tables using the following definition of an 
effective potential Ves, 



exp 



KfT(0) 



E exp 



Lt-Cs 



(80) 



We notice that this summation can be evaluated numerically to any desired accuracy. 

We calculated the thermodynamic averages after we had made visits to the array 
updating the phases and M visits updating the m's. Typically, if is small we used 
A^ = 4 and M = 1. In the opposite limit we used A^ = 1 and M = 8,10,.... This 
is so because our local updating algorithms for the m's have serious decorrelation time 
problems, due to the long range interaction among the charges. We typically found that in 
order to get reasonably small statistical errors, we needed to perform, in most cases, about 
Ameas = = 4096 measurements of the thermodynamical quantities, other times we took 
up to Ameas = = 8192 measurements. 

Once we have a long stationary string of values for the measured operators we calculated 



their mean values and uncertainties. We also have used the algorithm proposed in Ref. |^ 
for the efficient calculation of the helicity modulus. This method has a bias problem due 
to the last term in Eq. (|7TD . However, in the zero magnetic field case this problem is not 
present, since this term is identically zero. 



C. Results for f=0 

In this subsection we present the bulk of our Monte Carlo results. We have mostly 
calculated the helicity modulus in the small region and the inverse dielectric constant in 
the large regime, and both quantities in the intermediate region. 
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Most of the calculations we performed were for parameter values close to or at the 
experimental ones. In particular, the ratio between the self and mutual capacitances was 
kept fixed between the values Cg/Cm ~ 0.01 and 0.03, with the bulk of the calculations 
carried out for 0.01. We found that for the helicity modulus both values gave essentially 
the same results. Almost all of the calculations were done by lowering the temperature, in 



order to reduce the possibility for the system to be trapped in metastable states ||25[| . 

We have a clear physical understanding of the behavior of the system in the very small 
am limit, since this limit is close to the classical 2-D XY model. Moreover, we have the 
semiclassical calculation results, mentioned before, up to first order in am, which, as we 
shall see, agree very well with the Monte Carlo results. In this limit the results are solid 
because the discrete imaginary time path integral calculations converge very rapidly to the 
infinite limit. Therefore in this section we will discuss our numerical results for increasing 
values of a^. This will allow us to go from a well understood physical and calculational 
picture to the nonperturbative region of parameter space which is less understood. Here is 
where we will explore the limits of our numerical calculational schemes. The end result will 
be that a significant portion of the phase diagram can be understood. However, some of the 
most interesting intermediate regimes of the phase diagram are still very difficult to fully 
understand with our present calculational techniques. 

In Fig. ^ we show a typical curve for the helicity modulus as a function of temperature, 
in the small am limit. As a^ increases T flattens in the superconducting region. In order 
to calculate the transition temperature we used the fact that the critical temperature and 
the helicity modulus still satisfy the universal relation, 

T(T,) = - T,. (81) 
vr 

Based on the first order results from the semiclassical approximation analysis we know that 
this universal result is independent of am- In other words, we can determine the critical 
temperature by the intercept of T(T) with the line {2/n)T, as shown in Fig. |^. 

As can be seen from Fig. |^, at high temperatures and small am the asymptotic limit 
Lt- ^ oo is already reached for small L^. From Eq. ([T^), we can see that the parameter 
that determines this rate of convergence is 

(82) 



il3Ej)ara' 

The deep quantum limit is reached for a relatively large P ^ 1. This progression is shown 
in Fig. ^ where we plot the helicity modulus as a function of temperature for a relatively 
large am = 1.25, = Ly = 20, and three values of L^-. It can be seen that convergence is 
reached for P > 5, as found before in the self capacitive model in Ref. |]14| . 



As shown in Fig. |^, for a larger am the behavior changes and the departure from the 
Lr —>■ oo limit is manifested as a small dip in T at low temperature |l^ . As the temperature 
is lowered T shows an upward behavior. This is also seen in Fig. ^ from more extensive 
calculation for a still larger a^'s. This finite behavior can be understood in terms of a 
plane decoupling along the imaginary time direction. A way to see this is to notice that if we 
take both contributions to the action given in Eq.(^) as independent, both of them would 
yield a low temperature transition. If all the Lr planes are considered decoupled, then the 
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Josephson coupling would be [3/ Lr- Therefore in this case the N-S transition would happen 



at Tn_s ~ l/.t'T- On the other hand, if we only consider the second term in Eq. (15) we 
see that plane decoupling would take place at Tdeci ~ (ctrnZ-^r)- Now, if Tn_s > T^ecx which 
implies < 1, then only the first transition would take place. If however is large 
enough, we could have Tn-s < Tdeci, producing the observed dip in the hehcity modulus. 
This is seen in Fig. ^ where the helicity modulus is shown together with the inverse dielectric 
constant. There it can be seen that the dip starts at about the same temperature where 
becomes finite, signaling that the nonzero winding numbers have become relevant. In this 
region the fluctuations between planes have a small energy cost in the action given by Eq. 



As we increase am further we arrive at a point where at low temperatures, and fixed Lr, 
T goes to zero as shown in Fig. |^. To understand the nature of the low temperature phase 
we have computed the equal space imaginary time correlation function 

a(r) = ( COS (0(f, r) - 0(f, 0)) ). (83) 

We evaluated this function at three different temperatures, for = 1.75 and Lj. = 32. 
The results are shown in Fig. ^ where we see that the appropriate value for needs to be 
increased as T is lowered. More extreme are the two temperature results for = 2.5 and 



Lr = 64, shown in Fig. [1^. The upper curve corresponds to ksT / Ej = 0.36. As seen in Fig. 
H at this temperature a value of Lr = 64 is enough to reliably calculate the helicity modulus, 
for in this case the planes are correlated with a short decorrelation time. The lower curve 
has ksT/Ej = 0.1. At this temperature the helicity modulus is zero and the correlation 
function has a very short decorrelation time. These results show that the low temperature 
discontinuity is related to a decoupling of the planes along the imaginary time axis. 

This plane decoupling does not show any dependence on {L^, Ly) for the cases considered, 
and the curves shown in Fig. ^ are reproducible within the statistical errors for other values 
of Lx and Ly. Again we point out that upon increasing Lr the decoupling temperature 
moves closer to zero temperature. We should note that from the WKB analysis there is 
a critical value for a above which the superconducting state is no longer stable. So as we 
consider larger values of am, the superconducting state will become less and less stable. As 
we mentioned the simulations were performed while lowering the temperature. In contrast, 
in Fig.|TT] we show results from lowering the temperature for a^ = 2.75. In this case T 
reaches a zero for T < 0.2. We reversed the process increasing the temperature. Up to the 
last temperature calculated the results are consistent with having zero T for am = 2.75. The 
low temperature state arises from a decoupling transition between the imaginary time planes 
which leads to an ensemble of decorrelated planes. Maybe the planes could get recoupled at 
higher temperatures but, as already mentioned, our local algorithm would take too long to 
realign these planes so as to produce a coherent state. 

As shown in Fig. ^ the temperature where T has a sharp drop changes with the size of 
the system. To see if in the limit L,- — > oo we still have a finite low temperature transition, 
we tried to extract it from the data for three different a^ values by plotting them against 



1/ Lr as shown in Figs. and |T^. From these figures it appears that for am = 2.0 and 
2.25 there is a nonzero transition temperature in the oo limit. We used a jackknife 

calculation to estimate the infinite Lr temperature and we found 

Tdeci(am = 2, ^ oo) = (0.0183 ± 0.009)(Ej/A;b), 
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and 

Tdeci(am = 2.25, ^ oo) = (0.0067 ± 0.0025)(Ej/A;ij). 



The same type of calculation was done for am = 2.5. The results are shown in Fig. 0. Here 
we found that 

Tdcci(am = 2.5, L^^oo) = (-0.013 ± 0.005) {Ej/kn). 

Therefore from these estimates we surmise that the critical value for am ~ 2.5, which is 
larger than the one estimated from the WKB-RG analysis. However, our QMC calculations, 
in particular at low temperatures, are not precise enough to make a definitive determination 
of the critical am- 

The calculations of the S-N transition line seem to indicate that the superconductor to 
insulator zero temperature transition occurs at Om ~ 3, which is quantitatively different 
from the T=0 estimate |p. 

The evaluation of the inverse dielectric constant is considerably more complicated since 
the insulating region we have am > 3, needing larger values of at low temperatures. 
Moreover there are serious critical slowing down problems due to the long range charge 
interactions which worsen as the size of the system increases. 

We should point out that, in comparing with the purely classical case the quantum 
2-D Coulomb gas studied here has the extra complication of the n — cf) coupling, as seen in 
Eq. (|12[). This introduces an imaginary component to the action. Therefore we are forced 
to integrate the n's introducing the new variables {m} leading to the action given in Eq. 
(p!5|). This is the action that we used to perform the Monte Carlo calculations. 

We use the expression given in Eq. (^), which for finite k and = Ly gives \k\ = 27r/L^, 



as an upper bound for the inverse dielectric constant The technical problems mentioned 
above made the calculation of the dielectric constant less reliable than that of T. The results 
obtained from the Monte Carlo runs were too noisy to give us a quantitative estimate of 
the conductor to insulator transition temperature if there was one. Our quoted results give 
only tentative values for the transition line. 

The results are shown in Fig. |l|, where we also plotted the results of the Monte Carlo 
calculation of the normal to superconductor transition temperature as well as the experi- 
mental results from the Delft group [|5[. We have fitted a straight line to the first seven 
points in this line and used a jackknife calculation of Tc(am) for small am- We obtained 

= (0.9430 ± 0.0042) - (0.1800 ± 0.0040)am + O (al^ . (84) 
Ej ^ ^ 

The value of the slope is in good agreement with the semiclassical approximation result given 
in Eq. (^Tf ). The dashed line gives am = 2.8 at T=0 and joins the last QMC point to T = 0. 
The line is only a guide to the eye. We have not performed detailed calculations around 
am = 3 since the required values of makes reliable calculations too computationally 
intensive to be carried out with current algorithms and computer capabilities. 



D. Results for f=l/2. 

We also have performed a few calculations of the helicity modulus for the fully frustrated 
case f=l/2. The results of these calculations are shown in Fig. |l]. The experimental results 
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of the Delft group show a transition temperature for ctm = at (kBT/Ej) ^ 0.3 0, while 
the classical fully frustrated case has a critical temperature close to {ksT / Ej) ^ 0.5. Taking 
this into account we have rescaled the Monte Carlo results so that the calculated value for 
the transition temperature for = 0, f=l/2; coincides with the experimental result. 

Taking the experimental value of the critical temperature for the / = case, we per- 
formed a least square fit for the five smallest am's considered and found, using a jackknife 
calculation, 

(%^) ^ (0.9787 ±0.0070) - (0.256 ±0.017)a„, + O(a^). (85) 
V / /=o 

In the same way, but now for the / = 1/2 case, we found 

(0.3188 ± 0.0015) - (0.2929 ± 0.0066)0^ + 0(a^J. (86) 




These rough calculations show that the slopes of both curves are very close. This result 
can be compared with Eq. d^), which confirms that the first order correction in am to the 
critical temperature should not depend on the value of the magnetic field. A similar result 
for the equality of slopes was obtained in Ref. |T^, using a Monte Carlo calculation for the 
self-capacitive model. 



VI. SELF-CONSISTENT HARMONIC APPROXIMATION 

We need an alternative analytic approach, in principle exact for fixed and finite L^, to 
further understand the QMC results at low temperatures. This is important because of the 
strong Lt- dependence in the study of the QUIT temperature, and the analytic WKB results 
are only strictly valid at high temperatures and L^- = oo. In this section we use a variational 
principle to evaluate the free energy for the JJA within a self-consistent harmonic approxi- 
mation (SCHA). This approximation gives increasingly better results as the temperature is 
lowered. Previous SCHA calculations |^,|TD|,|TP[ did not explicitly include the charge degrees 



of freedom, which are of significant importance in the analysis presented here. As a bonus, 
we note that the SCHA developed here could be used as the basis for developing an alter- 
native QMC algorithm to study the model at low temperatures and intermediate a values, 
were both the WKB and the standard QMC analyses have problems. 

We start with the following decomposition of the Hamiltonian given in Eq. (^) 

H = Ho + {Hj- Hh}. (87) 

where Hq = {H + {Hh - Hj)}, and 

Hj = ^'^E E [l-cos(0(r,ri)-0(r,f2))], (88) 

''' T=l <ri,r2> 

Hh = ^'y: E [<l>{r,n)-<l>{T,r,)]\ (89) 

T=l <ri,r2> 
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In other words, we replace the Josephson Hamiltonian by a spin wave term and intro- 
duce its stiffness F as the variational parameter. F is of course the helicity modulus, 
but to emphasize that it is evaluated within the SCHA we us the F notation instead 
of T. The spin wave approximation does not have a phase transition at any tempera- 
ture so the variational calculation is going to give a vanishing F as the signature of the 
JJA transition. The variational free energy is given by Fy = Fh + {Hj — Hh) where 
I3Fh = —liaZn, and Zh = Tr {exp{— /Si^o}}, with the average <, > defined as usual by 
< A >H= Tr {Aexpj— /Si^o}} /Zh- We use the change of variables 



L-r-l 



/ 2 " / Tin 

^(r,f) = V^(0,f) + y— Xn(^) sin r 



r = l,2,...,L.-l, 



(90) 
(91) 



and perform the integrals over the variables Xn{^)- The result is 



ZH = Vdet[C]n E ^^Pj-^"^" 



xn 

r 



f m{r)=—oo 
2^ 0(0, f) 



Lr-1 

n 

n=l 



2LI [l - cos (f^) 



-LxLy /2 



X 



/o-^^^P^ ~o E 0(O,ri)N(fi,f2)0(O,f2) + Ej(O0(O,r) 



ri,r2 



(92) 



where 



h 



(2vr) 



'a 

f3r 



+ 



1 



1 



{[3EjT)gif3qJEjT/a 



X 



X 



E m(ri)0(fi, r2)m(f2). 



(93) 



ri,r2 



To obtain this equation we have taken Cg = 0, while the general case can also be treated as 
well, but since Cg <^ Cm in the experiment this assumption simplifies the calculations. We 
also have used the following definitions to write the previous equations. 



E ~ = E 0(n)O(fi,f2)0(f2), 



<Ti,r2> 



c 

N 



CinO, 



1-/1 l3qJEjV/C, 



o, 



(94) 

(95) 
(96) 



with O the lattice Laplacian operator, and g{*) and /(*) are functions defined in terms of 
a Matsubara sum and given in Appendix B. 

The details of the variational calculation of the free energy are presented in Appendix B. 
Here we discuss the main conclusions from the calculation. The results for am = 0, and 1 are 
shown in Fig. |l^. There we see that the am = case has the right low temperature linear 
T dependence. The result for am = 1-0 has an essentially flat low temperature behavior for 
F, which is due to quantum phase slips tunneling processes. In Appendix B it is shown that 
the helicity modulus can be expressed as 
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2 

with the fluctuations in the phase given by. 



T = exp\-U{A<pf)\. (97) 



,,,,2 sinhtejv/8^) 

^ ^ 2V r lcoshf/3Ejv/8^) -1 ' ^ ^ 



The classical limit corresponds to setting 0, which gives < (A0)^ >h{cI)= j^WJt- 

Using this result and Eq. (p^), at low temperatures T is given by F 1 — ksT /AEj + 0{T'^). 
This is precisely the same low temperature behavior obtained from a spin wave analysis 
in two dimensions [^. From < (A</))^ >h{cI) "we find the transition temperature within 



this approximation, Fc = ^ (^fy"), = 1/e, and kBT^^^Ej = 4/e ^ 1.472, which is an 
overestimate, (since the dimensionless 2-D critical temperature is T^"^ ^ 0.9). The problem 
with this approximation is extending the integration intervals from [— vr, tt] to (— oo, oo). As 
we have discussed this is a good approximation at low temperatures but it breaks down near 
the transition point. Our conclusions from this analysis are that for = Othe classical 
SCHA gives good results for low temperatures while it overestimates the critical temperature 
at higher ones. 

In the am 7^ case the quantization of the spin wave excitations leads to a non-vanishing 
result for < (A0)^ >h, given in Eq. (|9|). For (kBT/Ej) < 1 we get 

< (A0)2 (^1 + 2exp j-Z^Ej^S^I) , (99) 

j-ZJEj^to^l j , (100) 



To 



exp 



where Fq is the helicity modulus at zero temperature and it is the self-consistent solution 
to the equation Fq = exp| — y^^^j, with Fq ~ 1 — for ^ 1- The solution to 



this equation is shown as a function of in Fig. where we also show some Monte 
Carlo simulation results. The Fq result also presents a transition to a zero F state at 
a^(T = 0) = 32/e^ ~ 4.33, with a jump from Fg = 1/e^ to zero. Again the result of the 
SCHA overestimates the stability of the superconducting state. Both the extension of the 
integration intervals and ignoring the m's in these calculations are probably responsible for 
the deviations at large am- An interesting observation is that the result for Fc is exact up to 
first order in a^n- This is equal to the result we obtained from the WKB-RG analysis. Also 
surprising is that the first order correction to the critical temperature agrees with Eq. ( ^T]) 

'— am + 0(aj, (101) 



Ej / \ E.J \3z 



where z is the coordination number of the lattice. In Fig. |Tj we show F for a^a = 1.0 as 
a function of the temperature for increasing values of L^-, which should be compared with 
Fig. 1^. This result strongly suggests that the upward tendency of the helicity modulus at 
low temperatures seen in the QMC results may be an artifact of the finite nature of the 
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calculations. The origin of this increase is in the low temperature result for the finite Lr 
calculation of Eq. ( [B16|) 

<,.,).>„._|_(i__|^,o(r=)), ,102) 

SO that for finite Lr and at low temperatures the helicity modulus is given by 




E 



J 



+ 0{T'). (103) 



We have been able to explain the shape of the helicity modulus curves for low temper- 
atures, but Figs. ^ and ^ show that if > aj^, where ajj^ G (1.25, 1.75), then the hehcity 
modulus has a dip before it may go to one at low temperatures. So far we have ignored the 
contribution of the Discrete Gaussian Model (DGM) to the variational free energy which 
is a good approximation only for small am- We calculated the helicity modulus for the 
model ignoring the DGM and then including it as a continuous Gaussian model for finite 
Lt-, which is a good approximation for the effective coupling J^s <^ 1. For large the 
crossover point, which is when Jeff becomes soft, is seen as a finite dip in the helicity mod- 
ulus. Unfortunately, the T* found in the Monte Carlo calculations is much larger than the 
one given by Eq. (p7|) . It is apparent that the effective coupling Jeg does not contain all 
the contributions to the renormalization of the DGM due to the integration over the phases. 
To illustrate the nature of this crossover we performed several calculations of the helicity 
modulus for Lj. = 10, 20, and 40 with = 1. The results are shown in Fig. [T^. There we 
took ^^""^ = 6 for the crossover temperature. These results can be compared with those of 
Fig. |. 

The Monte Carlo calculations show that this dip occurs simultaneously with the rise 
in the inverse dielectric constant. This is a signal that the non-zero effective constant 
for the winding numbers becomes soft, making their contribution to the helicity modulus 
non vanishing. We note that the effect of a finite lattice, necessary for the Monte Carlo 
calculation, increases the softening temperature of Jes- On the other hand, the variational 
calculation for Cg = does not depend on the size of the lattice, therefore it can not capture 
these finite space size effects. 



VII. CONCLUSIONS 

We have presented a thorough study of the vs. T phase diagram for an array of 
ultrasmall Josephson junctions using a series of theoretical tools. One of our main goals 
was to perform these calculations for these arrays using experimentally realistic parameters. 
The model we used for the JJA is defined by a Hamiltonian that has two contributions, a 
Josephson coupling and an electrostatic interaction between the superconducting islands. 
The ratio of these two contributions was defined as am = (charging energy) / (Josephson 
energy). This was the important quantum parameter in our analysis. 

For convenience of calculation we derived different path integral formulations of the quan- 
tum partition function of the JJA. In the small limit we used a WKB-RG approximation 
to find the first order correction in to the classical partition function. The result of this 
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calculation was an effective classical partition function of a 2-D XY model type, where the 
coupling constant is modified by the quantum fluctuations. We used the modified renor- 
malization group equations for the 2-D model to find the superconducting to normal phase 
boundary. We also found that up to first order in the correction to the transition tem- 
perature was independent of magnetic field. One interesting finding from this calculation 
was the possible existence of a low temperature instability QUIT of the superconducting 
state. We found evidence for the QUIT, but the evidence is at the border of validity of the 
calculational approaches used. To have a definite theoretical proof of the existence of the 
QUIT one needs to have better algorithms and/or improvements in computer power. The 
results presented here are, however, rather encouraging. Of course the ultimate test will be 
furnished by experiment, and there too incipient indications of a low temperature instability 
have also been reported in 0. 

In the large limit we used a perturbative expansion in l/a^ to find an effective 
partition function for a quantum a 2-D Coulomb gas. This model shows a renormalized 
conducting to insulating transition as the temperature is lowered. We did not find a low 
temperature QUIT instability in the large insulating phase. 

We also performed extensive nonperturbative quantum Monte Carlo calculations of the 
JJA model. We concentrated our analysis on the helicity modulus T. This quantity is 
directly related to the superfluid density in the array. Using T we determined the super- 
conducting to normal transition boundary. We found good agreement between the critical 
temperatures obtained by QMC calculation and the semiclassical approximation. We also 
carried out a low temperature 1/Lr extrapolation analysis of the Tqujt and found evidence 
for Tquit 7^ for relatively large values of am = 2.0 and 2.25 but Tquit = for am = 2.5. 
These calculations have, however, a strong Lj. dependence and our QMC algorithm is not 
precise enough to completely ascertain the nature of the low temperature phase. Nonethe- 
less, the results found here together with the scant emerging experimental evidence for a 
low temperature instability yields further support for the possible existence of a QUIT. 

We also presented some QMC calculations of the inverse dielectric constant of the 2-D 
quantum Coulomb gas, in order to find the conducting to insulating phase boundary. We 
found that the present Monte Carlo path integral implementation of our model, that includes 
phase and charge degrees of freedom, does not allow us to make reliable calculations of this 
quantity. Our results for this transition are only qualitative. Further technical improvements 
are needed in order to make solid quantitative statements about the N — I phase boundary. 

To use a QMC calculation to prove or disprove the presence of a low temperature insta- 
bility in the superconducting state is not an easy task. However, these type of calculations 
give us upper temperature limits for the instability region. As far as our calculations could 
determine, the results for the superfluid density as a function of temperature are in rather 
good agreement with experimental findings in JJA except for the incipient data on the 
reentrant transition in the nonperturbative region of a^. 

To further understand the QMC results at low temperatures and as a function of L^, we 
have also implemented a self-consistent harmonic approximation analysis of the model, that 
includes phases and charge freedoms. We were able to make successful qualitative, and in 
some instances even quantitative comparisons between both calculational approaches. One 
of the conclusions from these calculations is that the general trend of the QMC results for 
T can be traced to the discretization of the imaginary time axis. For small a^, the decrease 
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of the helicity modulus at low temperature is clearly due to this effect. 

Among the most significant aspects of the results presented in this paper is the quan- 
titative agreement between our different calculational approaches and the corresponding 
experimental results in the superconducting-normal phase boundary with essentially only 
the measured capacitances as adjustable parameters. The existence of the QUIT is also a 
significant result of this paper, for it had not been studied in a model including realistic 
capacitances in the model. 
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APPENDIX: A 



In this appendix we present the derivation of the renormalization group equations in the 
insulating to normal region of the phase diagram up to first order in a^"^- We begin by 
writing the effective partition function from Eqs. (0), and (|54 



1 rl^h 



h Jo 



Zeff({m}) = [exp 
The phase average <></, is defined by 
<A> 



dTHj{{(t){f) + 0/(r,f) + {2t: / l3h)m{f)T}) 



(Al) 



IZ^'^f^^ A({0(f) + 0Kr,r)})exp 



From Eq.(^) the Josephson energy is 

Hj{{cP}) = -Ej E cos(0(fi) 



(A2) 



(A3) 



Inserting this equation into Eq. ([A2|) , and using the expansion given in Eq. ([5l|), we first 
see that the integrations over the {0}'s eliminate all the odd powers in Ej. Moreover up to 
the third term in Eq. (|5TD we have integrals of the form 



n 



27r 



d(f){r) 
2n 



exp 



0(fi) - 0(f2) + 



(A4) 



where rl ^ r2 and rs ^ f^, because they are pairs of nearest neighbors. 

Using this equation and the fact that Eq. ( ^41) is an even function of the 0/'s we end up 
with 
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^ofr({"i}) ~ 1 + ( ^ ) J J exp [- G(fi,f2;7-,r') 

^ <fl,f2> 

+ o{e: 



X COS (^{2tx / f3h)[m{fi) — m(r2)]( 
Here we have defined the correlation function 



r — T 



(A5) 



f 1 

ZiA -. J —oo 



exp 



j{T,r) =i 5{t- Ti) - 5{t - T2) 



(A6) 
(A7) 



The easiest way to perform this integral is to go to Fourier space, and using Eq. (|50D we 
can write 



0/(r,f) 



2 °° 

IIV^n(^)sin(z/„r), 

P''' n=l 



TT 



From this equation, Eq. ( [Aq ) can be written as 



(A8) 
(A9) 



f 1 00 POO 1 r , 

G{n, rl; ri, rs) = - In — n n / d^^e-n^f^^"^ exp 

I ^0 „=i f' ''-'^ 



n=l 



CO 

Sfl^n] = 2 II II i^nirs) {(/iZ^n/g)^C(f3,f4)}^„(f4), 
n=l ^^,7^4 

Therefore the correlation function is 

G(ri,r2;ri,r2) = -— ^ J2 — Jn(r3)C"^(f3, f4)j„(f4), 

n=l r3,r'4 



-q'PY. drj^ rfr'j(r,r-'3){g(r,r')C-i(r-'3,f4)}j(T',r4), 



where we have used the result (r < r') 



r, r 



Using Eq. ( [A7| ) we finally find 

G(fi,f2;ri,r2)=2g2/3J^'"'"'' 



1 sinfz/nr) sinfi/nT' 



(AlO) 
(All) 

(A12) 
(A13) 

(A14) 
(A15) 



(3h 



In - T2\ 

(3h 



C-\m - C-\\n - r,\) . (A16) 
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Inserting this equation in Eq. and noting that here \ fi — = \d\ we have 



X ^ exp< 

<r\,r2> 



2q^P \t - t'\ 



1 



T — T 



OiEr 



X COS [^{27T / l3h)[m{fi) — m{f2)]{T — r')^ -r '^v-^j; 
At this point we can perform the following symmetrical change of variables 

1 (t-t') 



Xi - 
X2 



V2 ph 

1 (r + rQ 

V2 ph 



(A17) 

(A18) 
(A19) 



Now we see that the integral over X2 can be calculated explicitly leaving only xi to be 
integrated. We find then 



Zeff({m}) ^ 1 + 



2q'P 



/•1/2 [ 

<ri,r2> ^ ' 



X cos {27r[m{n) -m{f2)]xi) + 0{Ej). (A20) 



Which is the result used in Eq. 



APPENDIX: B 



Here we give the details of the derivation of the effective free energy in the variational 
calculation, from which we obtain the expressions for < > h which was used to obtain 
the helicity modulus in Section [V^. 

We start by defining the following effective quantities for the m's and the 0(0, r) inter- 
action terms 



h 



-ttPEjT 



l-^)-h[Pq^EjT/C^ 



^m(fi)0(fi,r). 



ri 



(27r) 



pq^ 6 V lJ\ L, 



iPEjT)giPqJEjT/a 



X 



X 



m(fi)0(fi, r2)m(f2). 



(Bl) 



(B2) 



The functions h{A) and g{A) are given in terms of Matsubara frequency summations 



MA) ^ I f^V 'f: 



sin^(7m/L^) 



2\LJ [A2 + 2L2(1 - cos(7rn/L,))] 



1 — C0S(7™) 

1 — cos{TTn/Lr) 



2 /Lr 



Lr Lr \ A 



sinh((L^ - 1)A) + sinh(A) 
sinh(L^A) 



(B3) 
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A 



2 \Lr 



2 Lr-l 



- [A2 + 2L2(1 - cos(7™/L,))] [1 - cos(7™/L,) 



1 



A 



sinh((L^ - 1)A) fU - 1 



smh(LT-A) 



Lr 



1 / A 



6 \U 



{2Lr - 1; 



A = In n 



1 / A 



2 \Lr 



A 



\ 



A 



> . 



(B4) 
(B5) 



To complete the calculation of Zh we need to make one more approximation because the 
integral can not be explicitly calculated due to the finite limits of integration. Therefore we 
extend the limits of integration to the entire real axis. We expect that this approximation 
is good for large F or small T, which is our region of interest. For small F it would give an 
overestimation of the transition temperature. After extending the limits of integration we 
find that Eq. |9^ becomes 



Cri 



X 



1 + 



r Lr 

n 

.n=l 



l-h[(3q^EjTlC, 

LxLy /2 



-LxLy /I 



X 



2L2[1 -cos(7™/L^)]_ 
where the effective action for the m's is given by 

Sm = T^JeS m(fi)0(fi,r2)m(r2), 



n E expj-^i^l 



(B6) 



r m{r)=—Qo 



h"" 2 



(2vr) 



ri,r2 



(3q 



(3E,jT {(1 - l/L,) - h [(3q^E,jV/C^)Y 
A{l-h [Pq^EjT/C^ 



(B7) 



and O is the lattice Laplacian operator. The free energy in the harmonic approximation is 
therefore 



(3Eh _ PEnciJeE) , 1 
- +2' 



LxLy 



LxLy 



In 



Cri 



+ In 



l-hif3qJEjT/a 



+ 



Lr-l 

+ Ein 

n=l 



1 + 



2L2[1 -cos(7rn/L^)] 



(B8) 



The function is related to the partition function of the Discrete Gaussian Model (DGM) 
defined by 



PEDG{J) = -ln 



n E 

r m(r)=— oo 



J 



[m(fi) - m{f2)f 



<ri,r2> 



(B9) 
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We will ignore this term for the moment, since the discrete nature of its excitations makes 
it very small. We need to calculate the average of the harmonic part of the Hamiltonian, 
i.e. 



< Hh >h 



EjTLxLy 



dV 



H 13 dV 



l-hiPqJEjT/Q 



+ f (3qJEjT/a 



(BIO) 



where d joins nearest neighbors, and /(A) is given by a summation over Matsubara frequen- 
cies resulting in. 



/(A) 



(Lr - 1) 



sinh((L^ - 1)A) 



[1 + (A/L,) 2/2] sinh (L,A) 



(Bll) 



On the other hand we know that if the Hamiltonian is Gaussian, then the calculation of 
the average over the Josephson Hamiltonian is 



< a J >H= 2EjLxLy 



1 - exp i -i ([0(r, f + d) - 0(r, r)f) 



(B12) 



Therefore, the relevant term in the variational free energy is given by the average of the 
square of the lattice derivative over the harmonic Hamiltonian. Finally, we can write the 
variational free energy 



f3Fv (3Fh 



LxLy 



LxLy 



+ 2(3Ej 



l-exp^-i((A0)2 



< (A0)2 >H=< [0(r,f -0(r,r)]2 >^ . 



(B13) 
(B14) 



with < (A0)2 >H given in Eq. (|Biq ). 

The variational condition requires minimizing this function with respect to F. Taking 
the derivative of this equation, and using Eq. ( |B10|) , we find that 



(B15) 



F = exp{-i((A0f)}. 



F is then the normalized coupling constant of the spin wave Hamiltonian, i.e. the normalized 
helicity modulus. Eq. (^) gives the condition to identify the phase transition as the point 
at which F 7^ becomes a solution to the variational equations. The complete problem is to 
solve the combination of Eqs. ( pi4|) and (|97|). We can write Eqs. (|B14| ) as a function of F 



< (A0)^ > 



H- 



1 

2(3EjT 



1 - pEjJ2a^T 



h' 


[(3Ej^8a^rj 




1 - h 


[pEj^/Sa^f) 



+ f{(3EjJSa^T 



(B16) 
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We now study the behavior of the hehcity modulus as a function of in the continuum 



hmit, that is in the oo hmit. In this hmit we get 



cosh(A) 




< (Ac/))' >H 



sinh(A) 

1 

2 sinh(A)' ^ 2' 

sinh (^(3Ejy/8aj^T 



cosh (^(3EJ^/8a^] - 1 



(B17) 
(B18) 
(B19) 

(B20) 



We used Eqs. ( P15|) and ( P2Q|) to find the hehcity modulus for = and 1. The results 
of these calculations are shown in Fig. |I5[ This figure shows several interesting features. 
First the helicity modulus for = has the right shape, i. e., a linear T behavior 
near T = 0, and a jump at the transition temperature. The transition temperature is, 
however, overestimated. Note that the Om = 1 helicity modulus has a fiat behavior at low 
temperatures, which must be contrasted with the ctm = result. 
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FIGURES 



FIG. 1. The phase diagram for unfrustrated (f=0) and fully frustrated (f=l/2) arrays. Squares 
denote the experimental results from Ref. [5]. The other symbols give the QMC results including 
their statistical errors. The dotted lines for f=0 is a naive extrapolation of the QMC results to zero 
temperature, while the dot-dashed curve is the result of the WKB-RG calculation. In the f=l/2 
case we have rescaled the MC results so that the MC curve would coincide with the experimental 
result at am(/ = 1/2) = O.i 

FIG. 2. Renormalization group flow diagram. The discontinuous line indicates the vortex 
pair density as a function of temperature. See text for a discussion of the analysis of this RG flow 
diagram. 

FIG. 3. The helicity modulus T vs. temperature T for a^n = 0.5, = Ly = 32, and two 
different values of Lr = 6, and 8. This figure shows that for this value of a^i the results have 
converged to the infinite Lr limit already for Lr = 6. The line is {2kBT / Ejtt). The intersection of 
this line with T(T) is used to determine the normal to superconducting transition temperature. 

FIG. 4. The helicity modulus vs. temperature for Om = 1.25, Lx = Ly = 20, and three 
different values of Lr = 8, 16, and 32. This figure shows that a convenient convergence criteria for 
the Lr ^ CO limit is P = Lr/ {(3Ejara) > 5. 

FIG. 5. The same as Fig. |^ with am = 1.75, L^ = Ly = 20, and Lr = 16, 32 and 64. Again 
the criteria for convergence P > 5 works in this figure. Note that the departure from the infinite 
Lr limit is in the form of a small reentrant type dip in T as the temperature is lowered. 

FIG. 6. The helicity modulus vs. temperature for am = 2.25, L^ = Ly = 20, and 
Lr = 32, 48, 64, and 96. As in Fig. 5, T has a dip as the temperature is lowered and shows 
a reentrant-like behavior at lower temperatures. 

FIG. 7. T and inverse dielectric constant vs temperature for a^ = 2.25, L^ = Ly = 20, 
and Lr = 32. As the temperature is lowered, the drop in the helicity modulus is correlated with 
the rise in the inverse dielectric constant, signaling the decoupling between imaginary time planes. 

FIG. 8. T vs temperature for a^ = 2.5, L^ = Ly = 20, and Lr = 32, 48, 64, 80, 96, and 128. 
At this value of a^n the drop in the helicity modulus at low temperatures yields essentially T = 0. 
This abrupt drop is probably due to having a finite Lr- Note that the decrease of the T as T is 
lowered happens even before the plane decoupling temperature. 



37 



FIG. 9. The imaginary time correlation function C(r) vs. r, for am — 1.75, Lx — Ly — 20, 
and L-r = 32 for T = 0.35, 0.30, and 0.18. As the temperature is lowered, the correlation time 
decreases until at (kBT/Ej) ~ 0.3 the correlation time is shorter than one lattice spacing in the 
time direction. At this point the value of L^- is too small to produce a convergent calculation of 
the helicity modulus. 



FIG. 10. The imaginary time correlation function C(r) vs r for Om = 2.5, = Ly = 20, and 
Lr = 64, at (khT/Ej) = 0.36 and 0.10. As in Fig. |^ this correlation function for temperatures 
below the decoupling temperature is consistent with having zero decorrelation time. 



FIG. 11. The helicity modulus vs. temperature for a^n = 2.75, Lx = Ly = 20, and Lj- = 64. 
The upper curve was obtained while lowering the temperature. The lower curve corresponds to 
increasing it. The disparity in the results is most likely due to the difficulty in establishing phase 
coherence in the 64 equal time planes as T is increased, above the plane decoupling temperature. 



FIG. 12. Estimate of the plane decoupling transition vs l/L^ extracted from results similar 
to those of Fig. | for = 2.0, Lx = Ly = 20 and L^ = 48, 64, 80, 96, 128. The line from a least 
squares fit suggests a non zero decoupling temperature in the limit 1/Lr = 0. This figure is also 
consistent with the condition P ~ 5. 



FIG. 13. The SBiiie cLS Fig. |l 2| . Da,ta extrcictGd from Fig. ^ for dm — 

2.25 and Lx = Ly = 20. 

The line from a least squares fit also seems to suggest a non zero decoupling temperature in the 
limit L,- — > oo. This figure is also consistent with the P ~ 5 condition. 



FIG. 14. Same as Fig. 0. In this case, extracted from Fig. ^, for am = 2.5 and Lx = Ly = 20. 
The least squares straight line fit seems to suggest a negative decoupling temperature in the limit 

1/Lr = 0. 



FIG. 15. The spin wave stiffness vs T (related to the normalized helicity modulus) obtained 
from the self-consistent-harmonic approximation for a^a = and Om = 1- 



FIG. 16. The T = helicity modulus at as a function of a^, obtained from the SCHA. 
Extrapolated T = T results from the Monte Carlo simulations are shown as diamonds joined by 
a line. From this figure it is clear that the SCHA gives reasonable results at low temperatures and 
for small values of am- 



FIG. 17. The normalized helicity modulus F for am — 1.0 function of T for increasing 
values of Lr = 5,10,20, and 100. The extrapolated plot for L,- = cxd is shown as a solid line. 
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FIG. 18. The normalized helicity modulus T for ctj-Q — 1.0 and Lr = 10, 20, and 40. The 
crossover temperature T* is calculated using [L-,- / {(3Ej)a^^i) = 6, so that as L-^ increases the 
crossover temperature decreases. Here the solid line correspond to Lj- = 10, the dashed line to 
Lj- = 20 and the dot-dashed line to L,- = 40. See text for discussion of these results. 
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